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ABSTRACT 

We investigate a method to incorporate signal models that allow an additional fre¬ 
quency harmonic in searches for gravitational waves from spinning neutron stars. We 
assume emission is given by the general triaxial non-aligned model of Jones, whose 
waveform under certain conditions reduces to that of a biaxial processing star, or 
a simple rigidly rotating triaxial aligned star. The triaxial non-aligned and biaxial 
models can produce emission at both the star’s rotation frequency (/) and 2/, whilst 
the latter only emits at 2/. We have studied parameter estimation for signal mod¬ 
els using both a set of physical source parameters^ and a set of waveform parameters 
that remove a degeneracy. We have assessed the signal detection efficiency, and used 
Bayesian model selection to investigate how well we can distinguish between the three 
models. We found that for signal-to-noise ratios (SNRs) > 6 there is no significant loss 
in efficiency if performing a search for a signal at / and 2/ when the source is only 
producing emission at 2/. However, for sources with emission at both / and 2/ signals 
could be missed by a search only at 2/. We also find that for a triaxial aligned source, 
the correct model is always favoured, but for a triaxial non-aligned source it can be 
hard to distinguish between the triaxial non-aligned model and the biaxial model, even 
at high SNR. Finally, we apply the method to a selection of known pulsars using data 
from the LIGO fifth science run. We give the first upper limits on gravitational wave 
amplitude at both / and 2/ and apply the model selection criteria on real data. 

Key words: gravitational waves - stars: neutron - pulsars: general - methods: data 
analysis - methods: statistical 


1 INTRODUCTION 


Several searches have been performed for gravitational waves 
from known pulsars in data from the LIGO, GEO600 and 
Virgo gravitational wave detectors (Abbott et al.|2005 2007 


2008 2010 Abadie et al.|20lT Aasi et al.|2014 ). These rely 

on the known phase evolution of the pulsars from electro¬ 
magnetic observations (e.g. [Manchester et al.|200^ ) to allow 
long duration (of order a year) coherent searches for sig¬ 
nals from them in gravitational wave data. Unfortunately 
no signal has yet been seen, but interesting upper limits 
on gravitational wave amplitude have been produced, and 
for two pulsars (the Crab and Vela pulsars) the “spin-down 
limit” has been beaten ( Abbott et al.|[2Q0^ [Abadie et al. 
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2011). One of the principal previous methods used for these 


searches (Dupuis V Woan 2005) has focussed on parame¬ 
ter estimation and the setting of upper limits, but has not 
provided any measure of detection confidence. 

Previous gravitational wave searches targeting known 
pulsars have assumed gravitational wave emission at a sin¬ 
gle frequency, taken to be twice (or very close to twice) 
the spin frequency. However, there are reasons to consider 
slightly more general waveforms. In this paper we consider 


the model proposed by Jones (2010), dealing with steadily 


rotating triaxial stars. Specifically, Jones (2010) considered 


a star containing a pinned superfluid. Such a star can ro¬ 
tate steadily about an axis that does not coincide with the 
principal axis of the solid crust, and will generically emit 
gravitational radiation at both the spin frequency / and at 
2/. We term this the triaxial non-aligned case. This con¬ 
trasts with the ‘standard’ scenario, considered in almost all 
gravitational wave searches to date, of rotation about a prin¬ 
cipal axis, which emits only at 2/. We term this the triaxial 
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aligned case, and it can be regarded as a special case of the 
triaxial non-aligned case. Another special case is that of a 
biaxial star, where two moments of inertia of the star are 
equal. The waveform in this case is identical to that of a bi¬ 
axial processing star of the sort considered by [Zimmermann] 


& Szedenits (1979), which also produces gravitational waves 


at two frequencies. However, precession generically results in 
a modulation in the electromagnetic signal produced by a 
pulsar (see e.g. Jones Andersson|| 2001 ), something that is 
not clearly observed in the pulsar population. In contrast, 
in the model of Jones (2010), there is emission at / and 2/ 
even in a steadily spinning star. The attraction of this model 
is that such emission, at both / and 2 /, might be being pro¬ 
duced by any of the known pulsars, without leaving any tell¬ 
tale signature in the radio pulsations. It is therefore clearly 
of interest to understand the issues that arise when carrying 
out gravitational wave searches for such double-component 
signals. 

In this paper we study how different parameterisations 
of the model affect the estimation of signal parameters and 
the astrophysical information that can be extracted. We also 
discuss applying Bayesian model selection to assess the de¬ 
tection of signals from these sources and perform compar¬ 
isons between the different signal models. A similar study 
has been performed by Bejger & Krolak (2014) although 
there analysis was based on a maximum likelihood approach 
to parameter estimation. We use the methods we have de¬ 
veloped to analyse data from LIGO’s fifth science run (S5), 
setting upper limits on the emission at both / and 2 / for 43 
known pulsars. 

The plan of this paper is as follows. In Sectionj^we give 
a brief description of the neutron star model and waveform, 
confining the details to Appendixj^ In Section|3|we describe 
the Bayesian methodology we employ. In Section]^we briefly 
look at the shape of the parameter probability distributions 
for two different signal parameterisations. In Section we 
show how these Bayesian methods can allow us to distin¬ 
guish between the three different sorts of signals described 
above. In Section we present the results from applying a 
search for gravitational wave emission at both / and 2 / in 
LIGO data. We summarise our findings in Section 


2 THE MODEL 


In this Section we describe the physical model and gravita¬ 
tional wave emission from our triaxial star. In Section [ tH we 


use the original parameterisation of Jones (2010), while in 
Section [ 2 ^ we use an alternative simpler set of parameters, 
as described in Jones (2015). The ranges of the relevant pa¬ 
rameters are given in Section|2.3| again based on the analysis 
of|Jones|([2015l. 


2.1 The signal written in terms of source 
parameters 


Here we recap the physical model given in Jones (2010). 
The neutron star is triaxial, with a moment of inertia ten¬ 
sor whose principal components are (/i,/ 2 ,/ 3 ). Because of 
superfluid pinning, it can rotate about an axis, fixed in the 
inertial frame, that does not coincide with any one of these 
principal axes. This gives rise to gravitational wave emission 


at both / and 2/. The signal in a detector at the rotation 
frequency (/) is d Jones|[ 2 ^ 

(t) t) sin^cos^|/ 2 i sin 2 Asin^cos 0 (t) + 

(/21 cos^A — / 31 ) sin2^sin0(t)| 

— Fx ('0, t) sin^|(/ 2 i cos^A — / 31 ) sin 2^cos0(t) — 

/21 sin2Asin^sin0(t)|, (I) 

and the signal at twice the rotation frequency ( 2 /) is 

=2T+(?/;,t)(I + cos^f.)| [/ 2 i(sin^A — cos ^ A cos ^^) — 

/31 sin cos + I 21 sin 2A cos 0 sin 20(t) | 

— 4Tx (' 0 , t) cos^|/ 2 i sin 2 Acos^cos 2 (/)(t) — 

[/21 (sin ^ A — cos ^ A cos ^0) — I 31 sin ‘^0] sin 2(j){t) |. 

( 2 ) 


The polarisation factors T+ and Fx depend upon the po¬ 
larisation angle 'ip of the source. They also depend on the 
position of the source on the sky. We have not explicitly la¬ 
belled this dependence as these parameters would be known 
for a targeted gravitational wave search. The angle l is the 
inclination angle of the star’s spin vector with respect to the 
observer. 

The evolution in phase is generated by the rotation 
of the star, so that = 27v f{t') dt' + 0 o where f{t) is 

the frequency evolution and 0 o the phase at to- In prac¬ 
tice, for targeted gravitational wave searches, f(t) will be a 
known function (known e.g. from radio pulsar observations), 
and so we will simply write 


0(t) = Qt + 00, 


( 3 ) 


treating Q = 27r/ as a constant. 

The constant angles {0, 0o, A) are the Euler angles that 
specify the orientation of the star with respect to the iner¬ 
tial frame (at time to)- Here we have used A to replace the 
‘0’ parameter in Jones (2015) to avoid confusion with the 


standard use of 0 for gravitational wave polarisation angle. 
The parameters I21 and I31 are measures of the asymmetry 
in the moment of inertia tensor, with factors of the angular 
spin frequency Q and distance r absorbed for convenience: 


J -21 = -, -/31 = -- (,4j 

T r 

Putting all of this together, and assuming that the sky po¬ 
sition and spin frequency are already known, we have a set 
of seven source parameters: 


^source "(^5 05 -^21 5 -^31 5 ^5 005 Aj-. (5) 


We term this general case the triaxial non-aligned model 
of a spinning neutron star. There are two special cases that 
we will single out. The first is a triaxial star spinning about 
a principal axis. This can be obtained from Eqns ^ and 
§ by setting ^ = 0 ; there is then emission only at 2 /. 
We term this the triaxial aligned case. The second special 
case is the biaxial case, where two of the principal com¬ 
ponents of the quadrupole moment tensor are equal. This 
can be found by setting I 21 — 0 , and produces emission at 
both / and 2/. Note that, physically, this is slightly differ¬ 
ent from the relatively well-known precessional motion of a 
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biaxial star (see e.g. Zimmermann Szedenits||1979 Jones 


Andersson||2002 ), as the latter has an additional slow ro¬ 


tation, superimposed about the symmetry axis. However, 
the time variation of the mass quadrupole, and therefore 
the corresponding gravitational waveforms, are identical in 
the two cases, so all of the discussion of the biaxial case in 
this paper applies also to the biaxial precession waveform. 
Nevertheless, it should be remembered that in the preces¬ 
sion case, there can be modulation in the observed electro¬ 
magnetic pulsation frequency, and the time average of this 
electromagnetic pulsation frequency can be offset from the 
gravitational wave frequency; see Jones & Andersson (2002) 
for a detailed discussion. 


Table 1. Parameter ranges for the waveform parameters for the 
three signal models. 


Models 


Triaxial aligned 

Biaxial 

Triaxial non-aligned 

C 21 

— 

/^max /^max 

^21 W 21 

0, C™ 

C 22 

0, C™ 

/^max /^max 
'-^22 5 ^^22 

0, C™ 

^21 ('■ads) 

— 

0,27r 

0,27r 

< 1>^2 (rads) 

0,27r 


0,27r 

0 (rads) 

0,7r/2 

0,n/2 

0,7r/2 

cos L 

-1,1 

-1,1 

-1,1 


2.2 The signal written in terms of waveform 
parameters 


As previously shown by one of us (Jones 2015| ), and also 
explained in Bejger & Krolak (2014), the physical source 
model, specified by the seven parameters of Eqn. , con¬ 
tains a degeneracy. If we instead express the model as com¬ 
plex harmonic amplitudes we hnd that Eqns. ^ and § can 
be rewritten as 


h^{t) = — ^E+(' 0 , t)C2i sin ^cos ^cos + ^21^ — 

^Fx{ip,t)C2i Sint sin (6) 

and 

= - E+('0, t)C 22 [l + cos^^] cos (2(j){t) + ^^ 2 ^- 
2Fx ('0, t)C 22 cos L sin (^20(t) + ^>^ 2 ) • (7) 


There now appear two amplitude-like parameters C 21 and 
C 22 with corresponding phase parameters 4>2i and 4>22- As¬ 
suming that the sky location and spin frequency are known, 
we can identify a set of 6 waveform parameters^ one fewer 
than in the case of the source parameters: 

^waveform = <^21, ^22,4>2l, 4>22} • (8) 


When expressed in terms of these waveform parameters, a 
problematic degeneracy is removed, as we will illustrate in 
Section ^ below. Comparing with the source parameters of 
Eqn. (15), we see that the two angles giving the ori¬ 

entation of the star’s spin axis relative to the observer are 
common to both sets. There is in fact a rather complicated 
algebraic relationship between the hve remaining source pa¬ 
rameters (/ 21 , / 31 , 00, A) and the four remaining waveform 

parameters (C 21 , C 22 , 4>2i, 4>22)- This relation is derived 


Jones ( 2015| , and reproduced in Appendix where we 


also summarise the form that the waveform parameterisa- 
tion takes when specialised to the triaxial aligned and bi¬ 
axial cases. As shown in 


Jones 


(2015), the parameters C 21 
and C 22 are basically the (moduli) of the (complex) mass 
quadrupole scalars that describe the quadrupolar compo¬ 
nent of the mass distribution of the rotating star, with a 
factor of ord er absorbed for simplicity. (Note that in 


Jones 


(2015) these quantities are denoted by C 2 i,C 22 ). 


Note that in this analysis we are assuming a search for 
gravitational wave signals from known pulsars, or sources 
where a signihcant gravitational wave signal has already 
been found. This means that, rather than using the wave¬ 
forms as written in Eqns ([Tlf^ , or Eqns (IM^ , we can remove 


the oscillations that take place at the relatively high frequen¬ 
cies / and 2/, using the heterodyne method of [Dupuis 


Woan (2005) to give instead a pair of narrow-band complex 


times series. We do this by using the known phase evolution 
of the signal 0(t), multiplying by for the /-band and 

^-i 2 cf){t) 2/-band. This heterodyning, and subsequent 

low-pass hltering, leaves a signal model for the / and 2/ 
streams of 




.C 2 


^/(^) —-sin ^cos (0, t) sin 

( 9 ) 

and 


h2f{t) — -^E+('0,t)[l+cos^0e'^22_^^(y'22Ex('0,t)cosm'^22 

( 10 ) 

when written in terms of the waveform parameters. A similar 
heterodyning can be applied to the waveform when written 
in terms of the source parameters (see |Gill|2012 ). 


2.3 Parameter ranges 


In order to carry out our analyses, we need to choose sen¬ 
sible ranges in both the source and waveform parameters, 
for each of the triaxial aligned, biaxial and triaxial non- 
aligned models. The choice of ranges in these parameters 
turns out to be rather subtle, and is described in detail in 


Jones (2015). Basically, the source parameterisation, and. 


to a lesser extent, the waveform parameterisation, contain 
various discrete degeneracies, where changes in some com¬ 
bination of angle and/or amplitude parameters leaves the 
detected waveform h{t) invariant. This allows the ranges in 
these parameters to be reduced as compared to one’s initial 
expectations, with there being several options as to how the 
parameter space is reduced. The choices give n in Tables 
and represent one of several possibilities (see Jones|20T5 
and are the ranges we have used for the subsequent analy¬ 
ses presented in this paper. Any point in these ranges can 
be mapped into another part of the full parameter range 
that gives an identical waveform through the transforma- 


f 


tions given in Jones (2015). This enables signal parameter 


estimation and evidence evaluation to be performed using 
this minimal range, but for posteriors to then be mapped 
into the full range, if so desired. 


3 BAYESIAN METHODOLOGY 

Eor this analysis we want to be able to compute probability 
distributions for source and waveform parameters, and also 






























































4 M. Pitkin et al 


Table 2. Parameter ranges for the source parameters for the 
three signal models. 


Models 


Triaxial aligned 

Biaxial 

Triaxial non-aligned 

hi 

— 

0,/3™ 


hi 

n rmax 

U, ±21 

— 

0, 7™ 

00 (rads) 

0, TT 

0,27r 

0,27r 

A (rads) 

— 

— 

0, TT 

cos 6 

1 

0,1 

0,1 

0 (rads) 

0,7r/2 

0, TT 

0,7r/2 

cos L 

-1,1 

-1,1 

-1,1 


to compare models (noise-only verses triaxial non-aligned 
verses biaxial verses triaxial aligned). In Bayesian method¬ 
ology the standard way to compute probability distributions 
for unknown parameters is make use of Bayes theorem for 
the posterior probability distribution 


p{e\d,M,i) 


p{d\9,M,I)p{9\M,I) 
pid\M, I) 


( 11 ) 


where p(d\6, M, I) is the likelihood of the data d given model 
M and background information /, with a set of parameters 
0, p{0\M, I) is the prior on the parameters, and p{d\M, I) is 
the evidence^ or marginal likelihood (in this paper we will use 
the term evidence throughout for consistency), of the data 
given the model. The evidence is the factor that normalises 
the posterior probability density. It is given by 

p(d|M,/)= f p(d\G,M,I)p(G\M,I)de. (12) 
Je 

To compare models, we can calculate the Bayes factor, 
or odds ratio^ between competing models. To this end, note 
that for any model we can calculate its posterior probability 
as 


p{M\d,I) 


p{d\M,I)p{M\I) 

pid\I) 


(13) 


It is hard (maybe impossible) to calculate the normalisation 
factor p{d\I) as you have to know all alternative models and 
marginalise over them, but we can still compare posterior 
probabilities between models provided they use the same 
data. We can compute the Bayes factor, or odds ratio O 
(which we will use from here onwards) between two models, 
as e.g. 


p{Mi\d,I) p{d\Mi,I)p{Mi\I) 
piM2\d,I) p{d\M2,I)piM2\I)' 


(14) 


Note that the normalising factor p{d\I) has canceled out. If 
there is no known prior preference between the two models 
then the ratio p{Mi\I)/p{M 2 \I), the ratio of the prior odds 
between each model, can be set equal to unity. In this case, 
the odds ratio is just the ratio between the evidences, given 
by Eqn. (12), of the two models. We will adopt this view¬ 
point here, so all odds ratios will be calculated as the ratio 
of evidences. 

In the analyses performed in Sections and the like¬ 
lihood function p(d|0,M,/) we use is the Student’s t likeli¬ 
hood given in [Dupuis Woan (2005). This likelihood as¬ 
sumes that the noise in the data is stationary (over the 
dehned length of time) and Gaussian, but with an un¬ 
known noise standard deviation that has been analytically 


marginalised out. However, for the analysis of real data in 
Section [ 6 ] we have instead estimated the noise level for each 
data point and therefore use a Gaussian likelihood function 
in that section. The reason for this difference in likelihood 
function is that for real data it is more efficient to produce 
our processed data set at a lower sample rate and with the 
noise already estimated, which makes the Gaussian likeli¬ 
hood more appropriate. However, for large numbers of data 
points the two likelihoods will be very similar. 


3.1 Priors 

To compute evidences and posterior probability distribu¬ 
tions we must also explicitly dehne our prior probability dis¬ 
tributions. For the azimuthal-type angular parameters, and 
uniform in the cosine of the polar-type angular parameters, 
the least informative prior is a uniform prior dehned within 
their allowed ranges. So, given the ranges in Tables and 
the prior on the angles in the waveform parameterisation, 
assuming the triaxial non-aligned model, are 


p(^ 21 , ^ 22 ,'ip, cos l\M, I) 


const, if 0 ^ ^21 ^ 27r 

and 0 ^ ^22 ^ 27r 
< and 0 ^ ^ 7 r /2 

and — 1 ^ cos ^ ^ 1 ; 
^ 0 otherwise. 


(15) 

whilst in the source parameterisation, assuming the triaxial 
non-aligned model, it is 


p(0o, A, cos cos ^|M,/) = < 


0 


"const, if 0 ^ 00 ^ 27r 
and 0 ^ A ^ TT 
and 0 ^ cos ^ ^ 1 
and 0^0^ 7 r /2 
and — 1 ^ cos ^ ^ 1; 
otherwise. 

(16) 

Equivalents of these priors for the required parameters are 
used in the triaxial aligned and biaxial cases. 

We will use priors on the amplitude parameters that are 
uniform within a range dehned by the limits in Tables and 
These limits on the priors vary for the different model 
types, as described in Jones (2015). For the waveform pa¬ 


rameterisation in the triaxial non-aligned case our prior is 


^(<^ 21 , 0221 ^^,/) = 


rconst, if 0 ^ C 21 ^ C 2 T" 
and 0 ^ C22 ^ C^2 
0 otherwise. 


but for the biaxial case it is 


(17) 


f const. 


P(C21,C22|M,7) = 


and C 22 IC 21 ^ 0 ; 
otherwise. 


(18) 

For the source parameterisation in the triaxial non-aligned 
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case we choose to use a prior on the amplitudes given by 




{ const, if 0 ^ /31 ^ 

and 0 ^ /21 ^ /2T" 
and hi > hi] 

0 otherwise. 


(19) 

These priors on the amplitude parameters are uniform 
largely for convenience and simplicity rather than through 
a physical motivation. This is consistent with the uniform 
priors traditionally used in searches for gravitational waves 
from known pulsars, where uniform amplitude priors play 
a role of allowing a relatively high upper limit to be set 
by the likelihood, consistent with the data. However, when 
evaluating evidences a choice of uniform prior does have an 
influence, as doubling an allowed parameter range doubles 
the effective prior volume. However, we limit the effect of 
this in our analysis by assessing the distribution of odds ra¬ 
tios between signal and noise models empirically and basing 
thresholds on that empirical distribution. Additionally some 
of the influences of the size of prior volume cancel when com¬ 
paring signal models. It is worthwhile noting that evidence 
values produced using the minimal parameter ranges given 
in Tables and are equivalent to those that would be pro¬ 
duced using the full parameter space (or, e.g. just doubling 
the ranges). This is because the likelihood volume within 
the minimal range is exactly reproduced in each of the equiv¬ 
alent volumes within the total physical range, along with the 
prior volume increasing by the same factor. So, the increase 
in likelihood volume and prior volume cancel out. 


3.2 Nested sampling 


To calculate odds ratios we need to evaluate the evidence 


for each model, and Eqn. (12) shows this to involve multi¬ 


dimensional integrals. For some parameters, or likelihoods, 
the integral may be analytic, or for low numbers of dimen¬ 
sions it may be possible to evaluate it on a grid, but more 
generally, efficient numerical integration techniques must 
be applied. Here we use the nested sampling algorithm of 


Skilling (2006), in particular the implementation of it based 
on that developed by Veitch & Vecchio (2010) and available 


in the LALInference software library (Veitch et al. 2015). 


Nested sampling attempts to simplify Eqn. (12) into a one¬ 


dimensional integral that can be easily numerically calcu¬ 
lated. It samples a number of live points from the prior pa¬ 
rameter volume, calculates the likelihood at each point, finds 
the minimum likelihood Lmin point to add to the evidence 
integral, and then samples a new point with a higher likeli¬ 
hood than Lmin. This process is repeated until the integral 
is computed to sufficient accuracy. 

The analysis methods and models we have 
used have been incorporated into a code called 
lalapps_pulsar_paraineter_estimation_nested, which is 
freely available in the LALSuite software repositor}0 


^ http://WWW.Isc-group.phys.uwm.edu/daswg/projects/ 
lalsuite.html 


4 WAVEFORM VERSUS SOURCE 
PARAMETERS 


It is useful to look at some plots that illustrate the very dif¬ 
ferent nature of the waveform and source parameters. To do 
so, we can make use of the samples produced during nested 
sampling, by probabilistically drawing a subset that repre¬ 
sent the posterior probability distributions of the model pa¬ 
rameters, using either the waveform or source parameters. 
The distribution of samples for an individual parameter (or 
subset of parameters) represent the posterior probability for 
that parameter marginalised over all other parameters. This 
amounts to integrating Eqn. (11) over the prior ranges given 
in Tables and for the required parameter(s). In Figs. 
andj^the one-and-two dimensional posterior parameter dis¬ 
tributions are shown for the triaxial non-aligned model for 
an almost linearly polarised signal (cos ^ 0) with an SNR 

of 20, when recovered using the waveform and source pa¬ 
rameters respectively. Equivalent plots for an almost fully 
circularly polarised signal (| cos^| 1) are shown in Figs.|^ 
and^ We present results for these two extremes in inclina¬ 
tion angle to give the reader an idea of the range of different 
posterior probability distributions than can be obtained. 

From Figs. and it can be seen that the waveform 
parameters show a rather simple uni-modal probability dis¬ 
tribution. This is especially evident for the close-to-linearly- 
polarised signal, which shows the posteriors to be largely 
uncorrelated and Gaussian in appearance; as has been seen 
in previous triaxial aligned analyses, the extraction of pa¬ 
rameters for circular polarisations is slightly more difficult. 


due to increased correlations between the parameters (Pitkin 


2011). In contrast, the probability distributions in the source 


parameter space, shown in Figs. and show a large 
amount of structure (as originally observed in Gill 2012). 
As described in Appendix ^ the five source parameters 
(hi, hi,0,(j)o, X) can be related to the four waveform pa¬ 
rameters (C21, 4>2i, G22, ^22)- This leads to the source pa¬ 
rameters forming a highly degenerate and curved tube-like 
structure, rather than the much simpler form of the wave¬ 
form parameters. The full complexity of this tube-like struc¬ 
ture is probably being somewhat masked by our choice to 
only show its projected marginalisations in two dimensions. 


We note that non-negligible probabilities exist out to 
large values if hi ~ hi, due to degeneracies with other pa¬ 
rameters. For example, in Fig. hi and hi are truncated 
to show the bulk of the posterior, but would otherwise show 
a long tail along the diagonal of the hi versus hi joint 
posterior plot. This long tail in hi and hi is particularly 
prominent given our choice of a uniform prior in the am¬ 
plitude parameters. If we had used a prior uniform in the 
logarithm of the amplitude parameters then this tail would 
be greatly suppressed. 


Glearly, it will be much simpler to work with the wave¬ 
form parameters rather than source parameters. There is 
also the issue of computational speed. For a stochastic sam¬ 
pling technique such as nested sampling, the efficiency of 
the algorithm is greatly increased if new samples can be 
drawn from a distribution that closely matches the actual 
likelihood distribution. If the true distribution is smoothly 
varying, uni-modal, and relatively unstructured then it can 
generally be well approximated by a multivariate Gaussian. 
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<721 <722 (rad/7r) $^2 (rad/w) ^ (rad/7r) 


Figure 1. Marginalised posterior probability distribution plots of waveform parameters for an almost fully linearly polarised {cosl 0) 
signal with an SNR of 20 covering the minimal parameter ranges of Table ^ The cross-hairs show the true parameters of the simulated 
signal. All posterior plots have been produced with a modified version of the triangle .py python package ( [Foreman- Mackey et al.|2014| >. 


However, for more complex distributions such an approxi¬ 
mation becomes invalid. 

Indeed, comparisons in which the analysis code has been 
run on the same data, but using the waveform and source 
parameter spaces, show that to produce a similar number 
of posterior samples the former runs ~1.6 times faster than 
the latter in the case of no signal, and ~1.8 times fasteij^ 
for a signal with an SNR of ~20. This clearly shows the 

^ However, this speed difference can greatly increase when run- 


problems caused when sampling from likelihood functions 
with complex structure. 

We can therefore see that working with the waveform 
model rather than the source model is simpler, both in terms 
of the dimensionality of the parameter space and the shape 
of the likelihood function, and the required computations 
are faster in the waveform case. For these reasons, in the 

ning with larger numbers of live points to get better sampled 
posteriors. 
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Figure 2. Marginalised posterior probability distribution plots of source parameters for an almost fully linearly polarised (cos 6 0) 

signal with an SNR of 20 (the same signal as used in Fig.[^ covering the minimal parameter ranges of Table[^ The cross-hairs show the 
true parameters of the simulated signal. 


model comparisons that follow in Section we will work 
exclusively in terms of the waveform parameters. 


4.1 Parameter space mapping 

In the calculations described above, we have used ranges in 
the parameters that were as small as possible, i.e. we used 
the smallest possible sets such that one could be sure that 
if a triaxial non-aligned signal were present in the data, pa¬ 
rameters could be found that matched the signal; see |Jone^ 


(2015) for details. However, in the event of a detection, other 
parameters can be found that match the signal, in a way 
described using the transformations given in Jones (2015). 
Some of these other parameter sets correspond to physi¬ 
cally distinct stars. For instance, if one finds a signal with 
a polarisation angle '0, there will exist three other solutions 
with -0 values that differ by successive additions of 7r/2, and 
having different values for some other the other parameters; 
this degeneracy can only be broken if additional (probably 
electromagnetic) information is available. There will also be 
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Figure 3. Marginalised posterior probability distribution plots of waveform parameters for an almost fully circularly polarised (cos 1) 
signal with an SNR of 20 covering the minimal parameter ranges of Table ^ The cross-hairs show the true parameters of the simulated 
signal. 


other parameter sets that correspond to exactly the same 
physical star, and differ only in a trivial way, relating to 
how one chooses to label the three Cartesian axes that one 
lays down on the spinning triaxial body. It is instructive to 
fill out the full parameter space, to make clear that these 
degeneracies exist, and test the transformation rules given 
(12015). 


Jones 


In the case of the waveform model, we need only en¬ 
large the range covered by the polarisation angle -0, whose 
minimal range was '0 G [0,7r/2], and so covered one quarter 


of the full range -0 G [0, 27r] that So, to construct poste¬ 
riors in the full range we can randomly split the posterior 
samples into quarters, and each successive quarter can be 
mapped into the adjacent parameter volume by successive 
application of the transformations given in the Appendix of 


Jones (2015) 


'0 ^ '0 + 7r/2, 
+ TT, 

^22 ^22 ^■ 


( 20 ) 
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Figure 4. Marginalised posterior probability distribution plots of source parameters for an almost fully circularly polarised (cos 6 1) 

signal with an SNR of 20 (the same signal as used in Fig.[^ covering the minimal parameter ranges of Table[^ The cross-hairs show the 
true parameters of the simulated signal. 


For the source model the restricted range only covers a 
sixteenth of the full range (where cos^ G [—1,1], -0 G [0, 27r] 
and A G [0, 27r]). This means that the posterior samples have 
to be split between the sixteen volumes and more complex 
transforms used to map them as given in Jones (2015). If 
one wishes to map out this full parameter space, consid¬ 
erable care has to be taken when carrying out the trans¬ 
formations, particularly when selecting the correct roots of 
inverse trigonometric functions. For this reason, we give in 


Appendix [B] a outline of the procedure used here, written 
as a simple pseudo-code. 

The full posterior plots for the linearly polarised signal 
used in Figs. and based on this mapping of posterior 
samples are shown in Figs and In the waveform param- 
eterisation shown in Fig. ^ whilst the amplitude parame¬ 
ters and inclination can be unambiguously recovered, it is 
impossible through the gravitational wave signal alone to 
be able to distinguish the combination of initial phases and 
polarisation angle between the four distinct modes. In the 
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source parameterisation things become even more complex. 
Whilst /31 and cos l can be reasonably well determined the 
other parameters suffer from strong degeneracies. In partic¬ 
ular there are always combinations of parameters that allow 
I 21 to be close to zero or large, which means that it may 
only be possible to ever set upper limits on this parameter, 
even in the event of a detection. This implies that for the 
triaxial aligned model precise determination of the individ¬ 
ual physical parameters will not be possible even for high 
SNR sources. Only some highly correlated combination of 
parameters will be precisely determined. 

The complex degeneracies of the posteriors for the 
source parameters show that trying to estimate parame¬ 
ter uncertainties using the Gaussian approximation of the 
Fisher matrix would most likely lead to highly biased results 
even at high SNRs. However, the waveform parameterisation 
looks to be far more amenable to estimation using the Fisher 
matrix (as is done in Bejger Sz Kr61ak[|2014 ), provided that 
the minimal parameter space is used and the multi-modal 
degeneracies in the full parameter ranges are subsequently 
accounted for. 


5 MODEL SELECTION 


In this section we use simulations of signals and noise to eval¬ 
uate how well we can distinguish between noise-only data, 
and data containing a signal of the form of one of our three 
models (triaxial non-aligned, triaxial aligned, or biaxial). In 
all our simulations we have adopted a noise level for the 
/ and 2 / data streams based on the initial LIGO design 
sensitivity for the 4km LIGO Hanford Observatory (HI), 
i.e. the noise level is not equal between data streams and 
signal-to-noise ratios will be affected not just by the signal 
amplitude, but by their frequency. All simulations assume 
one day of heterodyned data sampled at a rate of one sam¬ 
ple per minute, as has been standard in previous searches 


(Dupuis & Woan 2005) and that the evidence evaluation 


uses a Student’s t likelihood function. For the advanced de¬ 
tectors the sensitivity curve shapes will be similar, but with 
the low-frequency edge being pushed to lower frequencies. 
Note that in all of our analyses, we assume that the signal 
phase evolution is tracked exactly, so that the frequency / 
is a known parameter. 

Also in this section we will often refer loosely to the 
natural logarithm of the odds ratio. In (9, as the ‘odds ratio’. 
This number has the convenience of being far more robust 
against issues of numerical precision when dealing with very 
large or small likelihood values. 


5.1 Noise-only simulations 


The odds ratio itself tells you how much one model is 
favoured over another, but computed odds ratios will be 
numerically different for different noise realisations. |Jeffreys| 
(1998) gives a rule-of-thumb table for assessing the signifi¬ 


cance of odds ratios, however this can also be approached 
empirically by using simulations to determine their distri¬ 
bution over realisations. In particular, by assessing the odds 
ratios empirically we can alleviate the effect of our large prior 
amplitude volume in shifting the distribution of odds ratios 
towards low values even when signals could potentially be 



Figure 7. The cumulative probability distribution of odds ratios 
comparing the three signal models to Gaussian noise for data 
containing only noise. 


seen. We will therefore follow this second path. The depen¬ 
dence of the odds ratios we calculate on the prior volume of 
the amplitude parameters probably shows that our choice of 
amplitude priors (uniform up to some maximum, and zero 
above the maximum) could be improved upon. One alter¬ 
native option is a prior that is uniform up to a lower maxi¬ 
mum value and then uniform in the logarithm of amplitude 
for larger values. However, the simple choice we have made 
has the advantage that the same prior can be used for all 
sources, no matter where they sit within the detector’s sen¬ 
sitivity curve, with the knowledge that the prior will easily 
bound the bulk of the likelihood in all cases. 

For our noise-only analysis we ran 2 000 simulations, 
with / drawn from a uniform distribution between 50 and 
700 Hz, and the / and 2/ data streams generated by drawing 
the real and imaginary components from Gaussian distribu¬ 
tions with zero mean and a variance set by the expected HI 
noise at the given frequenc}j^ For each simulation a random 
sky location was chosen, generated from a uniform distri¬ 
bution on the sky, as this determines the antenna pattern 
functions F+ and Tx of Eqns it and ( p^ . Note that, for the 
reasons outlined previously, when calculating odds ratios, 
we used signals written in terms of the waveform parame¬ 
ters rather than source parameters. For each simulation the 
odds ratios of evidence for each model versus that for Gaus¬ 
sian nois^have been calculated using nested sampling, and 
the cumulative probability distributions of all these values 
are shown in Fig.[^ Note that for the triaxial aligned model 
only the 2 / data stream was used. 

Fig.j^has three notable features. First, we see that the 


^ The variance is calculated as cr(t)^ = »S'n(/)/4At, with Sn{f) 
being the one-sided power spectral density at / and At being the 
time series time step, which is 60s for this analysis. 

For Gaussian noise the evidence is calculated by setting the 
signal amplitude to zero in the likelihood function. 
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Figure 5. Marginalised posterior probability distribution plots of waveform parameters for the same signal as used in Fig. but covering 
the full physical parameter ranges. The cross-hairs show the true parameters of the simulated signal. 


odds ratio for all three signal models is, for a large 

fraction of the simulations, small, confirming that noise-only 
data is normally found to be more consistent with noise than 
with a signal, as expected. Second, for a given cumulative 
probability, the odds ratios for the biaxial and triaxial non- 
aligned models are much smaller than for the triaxial aligned 
model. There are two related explanations for this: it is in¬ 
trinsically less likely for noise to conspire to imitate a signal 
with components at both / and 2/, as this would require the 
noise to produce signal-like disturbances at the two widely- 
separated frequencies, with correlated properties (e.g. with 


consistent values for the parameters l and '0); and, there 
will be an Occam factor at play, such that in the absence 
of any evidence for a particular signal the simpler triax¬ 
ial aligned model, which has fewer parameters and a corre¬ 
spondingly smaller prior volume than the biaxial and triaxial 
non-aligned models, will be favoured. In particular the large 
extra prior volume added by the extra amplitude parame¬ 
ter in the triaxial non-aligned and biaxial models play the 
predominant role in the shift between the odds ratio distri¬ 
butions. Third, the curves show that the biaxial model is 
slightly favoured over the triaxial non-aligned model. This 
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Figure 6. Marginalised posterior probability distribution plots of source parameters for the same signal as used in Fig. ^ but covering 
the full physical parameter ranges. The cross-hairs show the true parameters of the simulated signal. 


again can be accounted for through the Occam factor, due 
to the additional parameter, and therefore prior volume, re¬ 
quired for the triaxial non-aligned model over the biaxial 
model. 

From these distributions we can set an odds ratio 
threshold at which we favour one model over noise at a given 
false alarm probability. If we choose a false alarm proba¬ 
bility of 1% we find threshold odds ratios for each model 
versus Gaussian noise alone are —6.3, —13.0, and —13.5 for 
the triaxial aligned, biaxial and triaxial non-aligned cases 
respectively. We will use these thresholds for calculating de¬ 


tection efficiencies below. However, it is worth noting that 
analyses of data of different lengths of time, and/or com¬ 
bining additional detectors, would produce different values 
for the odds ratios. This means that a threshold needs to be 
calculated for a specific analysis and that the values above 
are only relevant if using one day of LIGO HI data sampled 
at a rate of once per minute. A different threshold would 
be required for the analysis of real LIGO data presented in 
Section!^ and indeed in Section [6.1| we demonstrate a dif¬ 
ferent, but related, assessment of detection significance. It is 
also worth noting that these simulations have used Gaussian 
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noise, whereas the distribution of odds ratios for real data 
would most likely be different. 


5.2 Signal simulations 

To assess detection efficiencies for signals described by the 
three different models we have generated simulations includ¬ 
ing these signals. In all three cases we drew the signal sky 
positions from a uniform distribution on the sky, and chose 
amplitude parameters to give a uniform SNR distribution 
between 0 and 20. We chose the distribution of angular pa¬ 
rameters to be uniform within the prior ranges given in Ta¬ 
ble In the triaxial non-aligned case we drew values from 
the source parameters, with I21 drawn from a uniform dis¬ 
tribution between zero and an upper range, whilst for each 
I21 value I31 was drawn from a uniform distribution between 
I21 and the same upper range, thus ensuring that /31 > / 21 . 
To obtain a uniform distribution in the overall SNR for the 
combined / and 2 / signal each of these pairs of parameter 
was re-scaled such that the final distribution was uniform in 
SNR between 0 and 20. 

In this analysis we define the SNR for a discretely- 
sampled complex signal, Xi, as 

( N p /2 

/5= - + ®(a;*)^] I (21) 

where a is the noise standard deviation (assumed constant) 
in both the real and imaginary components of the data. The 
coherent SNR for a combined / and 2/ signal is then just 
given by Pcoh = (p/ + 

We generated a set of 2 000 simulations including signals 
from the full triaxial non-aligned model, 2 000 simulations 
including signals from the triaxial aligned model, and 2 000 
simulations including signals from the biaxial model, with 
coherent SNRs between 0 and 20. To assess model compar¬ 
ison at much higher SNRs we generated further sets of 500 
injections of triaxial non-aligned signals with SNRs of 50, 
100 and 500. 

For each set of injections we calculated the odds ratios 
of signal vs. Gaussian noise for each of the three different 
models. We then used these to compute detection efficien¬ 
cies, and also to compute odds ratios between the three dif¬ 
ferent signal models. We now present the results for each 
type of injected signal. 


5.2.1 Detecting triaxial non-aligned signals 


Using the odds ratio threshold for the 1% false alarm prob¬ 
ability found in Section [5T] we can work out the efficiencies 
of each model for detecting an injected triaxial non-aligned 
signal. We can do this as a function of the total SNR, but 
also as a function of the SNR in both the / and 2/ streams 
individually. We can also use these odds ratios to compare 
the signal models against one another. 

Fig. [ 8 ^ shows the efficiency for detecting triaxial non- 
aligned signals as a function of the combined SNRs in the / 
and 2/ data streams (the shaded regions give 95% credible 
intervals calculated using the method of |Cameron|2011 ).The 
three curves correspond to using the (correct) triaxial non- 
aligned model, and using the (incorrect) triaxial aligned and 
biaxial models. The equivalent efficiencies as a function of 


the SNR in each individual stream are shown in Fig.[^ with 
Fig. 9a assuming the (correct) triaxial non-aligned model, 
and Figs. 9b andassuming the (incorrect) triaxial aligned 
and biaxial models, respectively. 

From Fig. it can be seen that a 95% detection ef¬ 
ficiency is achieved for the biaxial and triaxial non-aligned 
models for SNRs >5.5. The detection efficiency using the 
triaxial aligned model is systematically lower. The biaxial 
and triaxial non-aligned models recover signals with nearly 
equal efficiency showing that even though the injected triax¬ 
ial non-aligned signal contains an additional phase parame¬ 
ter, recovery with the biaxial model generally finds a nearly 
equivalent parameterisation to match it (see below). It is ob¬ 
vious that the triaxial aligned model will not detect signals 
for which there is very little power in the 2 / stream, which is 
the reason for its detection efficiency curve in Fig. [ 8 a| always 
being slightly lower than that of the triaxial non-aligned and 
biaxial models. This is seen as an obvious effect in Fig. |9b| 
However, signals with SNR > 6 in the 2/ stream are still 
well recovered. (In contrast, when assuming a triaxial non- 
aligned or biaxial signal, detection is possible providing the 
combined SNR is sufficiently large, as made clear in Figs. [9^ 
and respectively). This again shows that, even for sig¬ 
nals with the full triaxial non-aligned model parameterisa¬ 
tion, the 2 / component can still be well-matched with the 
triaxial aligned model. This good match is straightforwardly 
apparent when thinking in terms of the waveform parame¬ 
terisation as the signal in the 2 / component has exactly the 
same form for both models, whereas this equivalence is not 
so evident when using the source parameterisation. 

Using the odds ratio values we have calculated, we are 
also able to compare the Bayesian evidences between signal 
models rather than just comparing a model to noise. For 
the biaxial and triaxial non-aligned cases we can take the 
ratio of odds ratios that we have already calculated as the 
noise evidence terms will cancel out. However, for the tri¬ 
axial aligned case we need to include the evidence for there 
being no signal present in the / data stream as part of the 
signal model. Therefore the evidence for the triaxial aligned 
model becomes the product of the evidence for a signal at 
2/ and only noise at /. We take model 1 to be favoured over 
model 2 if 0\2 = 0\I02 > 1. For each model pairing (triax¬ 
ial aligned versus triaxial non-aligned, triaxial aligned versus 
biaxial and biaxial versus triaxial non-aligned) we have com¬ 
puted the percentage of our simulated signals for which the 
numerator model is favoured as a function of SNR, e.g. 50% 
means that on average either model is equally likely. Results 
for this can be seen in Fig. |10a| 

We see that up to SNRs of ~10 the simpler triaxial 
aligned model is more often favoured, whilst at greater SNR 
the triaxial non-aligned and biaxial models become more 
probable, i.e. the fact that the signal looks more like the tri¬ 
axial non-aligned and biaxial models starts to overcome the 
Occam factors disfavouring them. For signals with coherent 
SNRs of ~20 we see that the simpler triaxial aligned model 
is still favoured ~30% of the time. This is mainly a result 
of the distribution of our population of simulated signals for 
which, when the coherent SNR ~20, about a third of the 
signals have an SNR of < 5 in the / component. Therefore, 
in these cases the / component is providing very little addi¬ 
tional evidence for the triaxial non-aligned, or biaxial mod¬ 
els, and the simpler triaxial aligned model is being favoured. 
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Figure 10. The percentage of simulations favouring a particular model (see legends) when the simulation contains 0 the full triaxial 
non-aligned model, 0 the triaxial aligned model, and 0 the biaxial model, as a function of the signal coherent SNR. Note that in 0 
the (blue) triaxial aligned/triaxial non-aligned results are barely visible beneath the (red) triaxial aligned/biaxial results. 


A striking point to note is that out to SNRs of somewhere 
between 20 and 50 the biaxial model is favoured more than 
50% of the time over the true triaxial non-aligned model. 
As noted in Appendix |A2| the biaxial model is just a special 
case of the triaxial non-aligned model with the constraint 
that ^22 — 24>^i, so in some small fraction of the simula¬ 
tions this criterion will be fulfilled and the Occam factor will 
result in the biaxial model being favoured (see discussion 
below and Fig. 12). However, this does not account for the 
fraction of the cases that the biaxial model is favoured. To 
explain that we must look at another degeneracy: when the 
signal is circularly polarised ^22 and ^ become very highly 
correlated (see, e.g.. Fig. [^, and in these cases a combi¬ 
nation of ^22 and '0 can be found such that ^§2 ~ 24>^i. 
Therefore these will again look like biaxial signals and the 
Occam factor will start to favour them. This is demonstrated 
in Fig. m which shows the parameter posterior probability 
distributions for a triaxial non-aligned signal with SNR ~20 


and |cosf.|~0.5, even with |cosf.| some way from unity the 
correlation between ^22 and 'ip is still strong. This means 
that for SNRs ~20 about half the population will still be 
able to support the biaxial model. Indeed in the case of the 
signal in Fig. the biaxial model is very slightly favoured 
over the triaxial non-aligned model by a factor of ~ 1.3. Even 
at SNRs as high as 500 this effect still means that ~ 15-20% 
of triaxial non-aligned signal simulations favour the biaxial 
model. 


We can also see this effect by looking at three illustra¬ 
tive waveforms, each corresponding to signals with SNR of 
500, in the hrst case picking out a waveform where the biax¬ 
ial model is favoured by in the second case the biax¬ 

ial and triaxial non-aligned models are equally likely, and in 
the the third case the biaxial model is strongly disfavoured 
by a factor of The hrst case is shown in Fig. 12 


which happens to be a simulation in which ^22 ~ 24>^i and 
cos^ = 0.23, so the waveform is essentially a biaxial wave- 
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Figure 11. Marginalised posterior probability distribution plots of waveform parameters for an SNR ~20 signal with cos 6 = —0.49 
covering the minimal parameter ranges of Table ^ 


form and both the biaxial and triaxial non-aligned models 
recover the correct waveform almost perfectly - the param¬ 
eters are also recovered consistently for both models. In this 
case the triaxial non-aligned model provides no extra infor¬ 
mation, so the Occam factor means that the biaxial model 
is favoured. The second case is shown in Fig. |13| in which 
^22 ^ but cos^ = —0.93. In this case almost identical 

waveforms are recovered for both models, but due to the ^22 
and 'ijj degeneracy the recovered best fit parameters are not 
the same. The third case is shown in Fig. |14| where again 
^22 9^ 2$^i, but cosf. = 0.23. This shows that whereas the 


triaxial non-aligned model produces an excellent fit in both 
the / and 2/ data-streams the biaxial model sacrifices any 
attempt at a good fit in the 2/ data in favour of getting a 
very good fit in the higher SNR / data. 


5.2.2 Detecting triaxial aligned signals 

Using the triaxial aligned injections we have again calcu¬ 
lated odds ratios for each model. The efficiencies of detect¬ 
ing these signals for each of the models, using the 1% false 
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(a) 



(b) 



(c) 


Figure 8. Each figure shows the efficiencies for signal detection 
(with a 1% false alarm probability) based on the odds ratios recov¬ 
ered assuming the three different signal models. In the data 
contained triaxial non-aligned injections, in 0 triaxial aligned 
injections, and in B biaxial injections. 





(c) 


Figure 9. The efficiencies for signal detection (with a 1% false 
alarm probability) based on odds ratios for data containing tri¬ 
axial non-aligned injections, plotted as a function of the signal 
SNR in the individual / and 2/ data streams. The analysis in 
Q assumed the (correct) triaxial non-aligned model, while the 
analyses in 0 and 0 assumed (incorrectly) the triaxial aligned 
and biaxial models, respectively. 

















































































































































































































































































































































Dual-harmonic GWs from NSs 17 



Time (days) Time (days) 


2.0 

1.5 

1.0 

0.5 

/ 0.0 
-0.5 
- 1.0 
-1.5 


xl0~ 


Real 


Imaginary 


Data 

Tnaxial non-aligncd 

Biaxial 

Injection 




Figure 12. The real and imaginary waveforms for the / and 2/ 
data streams of an SNR 500 signal for which the biaxial model is 
favoured over the triaxial non-aligned model by The black 

points represent the simulated data, the overlapping red and blue 
lines show a distribution of waveforms drawn randomly from the 
posterior parameter distributions for the triaxial non-aligned and 
biaxial models respectively, and the dashed black line shows the 
injected waveform. In this case the lines are nearly impossible to 
distinguish on the plot, as a consistent waveform is recovered for 
both models. 


Figure 13. The real and imaginary waveforms for the / and 2/ 
data streams of an SNR 500 signal for which the biaxial model and 
triaxial non-aligned models are equally likely. The black points 
represent the simulated data, the overlapping red and blue lines 
show a distribution of waveforms drawn randomly from the pos¬ 
terior parameter distributions for the triaxial non-aligned and 
biaxial models respectively, and the dashed black line shows the 
injected waveform. In this case the lines are hard to distinguish 
between as a consistent waveform is recovered for both models, 
although in the / data stream there is a minor discrepancy. 


alarm probability thresholds given in Section [5.1| are shown 
in Fig. [Sl^ 

The triaxial aligned model is just a special case of the 
two other models and we see that, as it is the simplest model, 
the best efficiency is achieved when we just use that model. 
However, the efficiency increase is relatively small compared 
to more complex models, with those models giving almost 
100% efficiency for SNR greater than ~7. 

The main advantage of just performing the analysis as¬ 
suming the triaxial aligned signal model is speed. Using the 
waveform parameterisation the analysis in this case is ~2 
times faster than using the triaxial non-aligned model (the 
speed ratio does not vary significantly with SNR as the likeli¬ 
hoods in both cases are fairly simple). However, if the signal 
does contain significant power at / then, as we have seen 
in Figs [8^ and 1 9b I assuming a purely triaxial aligned signal 
model could lead to signals being missed. 

The signal model comparison is shown in Fig. |10b| This 
shows that when a purely triaxial aligned signal is present 
then that model is always favoured over the more complex 
models. As we have noted earlier the triaxial aligned model 
is just a special case of the other models, so this result is 
purely down to the Occam factor rather than it being a bet¬ 
ter fit to the data. We also see the Occam factor in play when 
comparing the biaxial and triaxial non-aligned models, with 
the biaxial model being favoured for the majority of signals. 
There is a slow trend towards the triaxial non-aligned model 
being more favoured. The reasons for this are not entirely 


clear, but a possible explanation is that the extra free pa¬ 
rameter in the triaxial non-aligned model will allow it to 
more easily accommodate the necessary lack of signal in the 
/ data-stream. 

5.2.3 Detecting biaxial signals 

Finally, using the biaxial simulations we have calculated 
odds ratios for each model. Fig. [8^ shows the efficiencies of 
detecting these signals for each of the models, using the 1% 
false alarm probability thresholds given in Section [5.1| The 
results are similar to those for the triaxial non-aligned injec¬ 
tions of Fig.|8a| with the triaxial aligned model being some¬ 
what less efficient than the biaxial or triaxial non-aligned 
models, due to its inability to detect signals with low SNR 
in just the 2/ component. Again we see that the efficiencies 
of the triaxial non-aligned and biaxial models are very sim¬ 
ilar, despite one of them (the latter, in this case) being the 
correct model. 

The signal model comparison is shown in Fig. |10c[ In 
common with the cases described above, at low SNR the tri¬ 
axial aligned model is favoured over the other two models. As 
the SNRs increase, the biaxial and triaxial non-aligned mod¬ 
els become more favoured over the triaxial aligned model, by 
approximately equal amounts, and are more often preferred 
for SNRs above about 10. The biaxial model is favoured 
over the triaxial non-aligned model by a constant amount of 
~90% over all SNRs except the very smallest ones. 















































































18 


M. Pitkin et al 



Figure 14. The real and imaginary waveforms for the / and 2/ 
data streams of an SNR 500 signal for which the biaxial model 
is disfavoured over the triaxial non-aligned model by a factor 
of The black points represent the simulated data, the 

overlapping red and blue lines show a distribution of waveforms 
drawn randomly from the posterior parameter distributions for 
the triaxial non-aligned and biaxial models respectively, and the 
dashed black line shows the injected waveform. In this case the 
biaxial model has to sacrifice any close fit in the 2/ data-stream 
over producing a very good fit to the higher SNR / data-stream. 


6 SEARCH IN REAL DATA 


We performed a search for gravitational waves from a selec¬ 
tion of isolated (i.e. non-binary system) known pulsars using 
data from the fifth LIGO science run (S5) ( [Abbott et al. 
2009). This is the first gravitational wave search targeted 
at known pulsars to include an explicit search for a compo¬ 
nent at the rotation frequency. We consider 43 pulsars that 
had previously been targeted using an analysis only sensi¬ 
tive to the triaxial aligned model with emission at twice the 
pulsars’ rotation frequencies ( Abbott et al.|[20T0 ). We used 
science mode data from the two LIGO Hanford Observatory 
detectors (HI and H2) and the LIGO Livingston Observa¬ 
tory detector (LI) covering 4 November 2005 to 1 October 
2007, which was processed into sets of discrete Fourier trans¬ 
forms using 30-min sections of data (called short Fourier 
transforms, or SFTs; as used in e.g. [Aasi et ^|2013| ). This 
gave a total of ~491 days of data for HI, ~497 days for 
H2 and ~392 days for LI. For each pulsar, we filtered the 
SFTs from each detector using a “spectral interpolation” 
routine ( Davies et al.|2015 ) to create two narrow-band com¬ 
plex data streams sampled at one per 30 minutes: one with 
the phase evolution at the pulsars’ rotation frequency re¬ 
moved, and the other with the phase evolution at twice the 
pulsars’ rotation frequency removed. For each data stream 
we also produced estimates of the time-varying standard de¬ 
viation of the corresponding noise. 

For each pulsar we performed parameter estimation and 
calculated the evidence for each of the triaxial aligned, bi¬ 


axial and triaxial non-aligned models. This allows us to per¬ 
form model comparison for each source, assess the detec¬ 
tion of signals and produce 95% credible region upper lim¬ 
its (bounded at zero) on the waveform model amplitudes 
C 21 and C 22 - We show the results, based on a joint anal¬ 
ysis of data from all three detectors, in Table As with 
the analysis in Abbott et al. (2010) we expect amplitude 
calibration uncertainties of 10% for HI and H2 and 13% 
for LI. The table includes results for the triaxial aligned 
model giving upper limits on the conventionally quoted grav¬ 
itational wave strain amplitude ho (which in this model 
is related to C 22 via ho = 2 C 22 ). These upper li mits are 
broadly consistent (within ~25^^ with those in 


Abbott 


et al. (2010), but are not identical due to differences in the 


processing pipeline (heterodyned, averaged and low-pass fil¬ 
tered ( Dupuis Woan]|2005| vs. Fourier-transformed) and 
the fact that here we have used a Gaussian likelihood for the 
data given our model (using a noise estimate for each sam¬ 
ple), rather than a Student’s t likelihood with its implied 
marginalisation over an unknown noise level. From Table 
we see that the triaxial aligned model is favoured by factors 
of ~e^^ in all cases. Given that the previous gravitational 
wave search reported in Abbott et al.| ( [20T0 ) found no evi¬ 
dence for triaxial aligned signals (purely based on examina¬ 
tion of posterior probability distributions of the estimated 
signal amplitudes), we can conclude that there is no evidence 
for biaxial or triaxial non-aligned signals either. However, 
below we will also assess the significance of the odds ratio 
for the triaxial non-aligned signal compared to noise, which 
appears in the results table. 

One advantage of using the waveform parameterisation 
over the source parameterisation is that C 21 and C 22 di¬ 
rectly represent the search sensitivity at / and 2/ respec¬ 
tively, whereas /31 and 1 21 contribute to both harmonics 
in a complicated way that cannot be disentangled. How¬ 
ever, it is important to note that we currently do not have 
a good understanding of the physical interpretation of C 21 
and C 22 . It is interesting to note from our results in Table 
that, although in many cases the detector sensitivities are 
better at / than 2/ (see e.g. Fig. 4 of Abbott et al.||201Q| ), 
we get smaller upper limits on C 22 for all bar two pulsars 
(J0024—7204G and J2322+2057). This because when ^ = 0 
the / signal is zero and C 21 can extend to arbitrarily large 
values, which creates a tail on the C 21 posterior probabil¬ 
ity distribution leading to the larger upper limits. The tail 
would be suppressed if a Jeffreys prior were used for the 
amplitudes, but would still be present at some level as the 
correlation is a feature of the waveform. Fig. shows an 
example of the posterior probability distributions for the 
triaxial non-aligned model for pulsar J0024—7204G. 


^ A notable outlier is J1748—2446ac for which our new result is a 
factor of 1.7 times better than that in [Abbott et al.jpOlOj ). This 
seems to be due to there being a wandering spectral line feature 
in the HI data close to 2/ for this pulsar, which the narrower 
bandwidth and noise estimation procedure of the spectral inter¬ 
polation method is able to veto, but that had artificially biased 
the noise level on the original result. 
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Figure 15. The posterior probability distributions for the waveform parameters of the triaxial non-aligned model found using S5 data 
for pulsar J0024—7204C. 


6.1 Results significance 

It is useful to have a way to assess whether the signal vs. 
noise odds ratio value for a particular pulsar is large enough 
to be considered a detection (or detection candidate). As 
mentioned earlier, the odds ratio itself tells us how much 
more probable the signal model is compared to noise given 
the data, but fluctuations in this value for different noise re¬ 
alisations and the effect of our large amplitude prior ranges 
mean that an understanding of the potential distribution 
of values is useful in making a detection decision. We did 
this in Section |5.1| using simulations of noise-only data to 


get a distribution of odds ratios when no signal was present, 
from which we could set a detection threshold, given a cho¬ 
sen false alarm rate. For data in which the noise is purely 
Gaussian this is straightforward, but with real data we need 
a way of producing the corresponding noise-only data with 
the same statistics to get a representative distributions of 
odds ratios. It is also useful to look at extra information 
such as the SNR of the maximum a posteriori recovered 
waveform. Assessing a particular search’s significance using 
an empirically estimated ‘background’ distribution of some 
detection statistic versus SNR is common for many searches 
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for transient gravitational waves, where in those cases the 
‘background’ is generated through many time-slides of the 
data. 

To estimate a noise-only distribution of log odds ra¬ 
tios for the triaxial non-aligned model versus noise for each 
pulsar we have made use of the same data as for the real 
analysis, but have ‘scrambled’ the data by randomly shuf¬ 
fling the time order. This preserves the same noise statistics, 
but would completely de-cohere any signal present, essen¬ 
tially giving us random realisations of the data. For each 
pulsar therefore, we shuffled the data 100 times and calcu¬ 
lated the log odds ratio for the triaxial non-aligned model 
versus Gaussian noise, whilst also recording the SNR of the 
maximum a posteriori recovered signal model. We calculated 
the correlation matrix of these 100 pairs of log odds ratios 
and SNR and then, assuming their distribution is a bivari¬ 
ate Gaussian, calculated a set of probability contours in the 
odds ratio-SNR plane for these. The actual odds ratio and 
SNR for the pulsar (obtained using the non-scrambled data) 
can then be placed on the plot, and its location relative to 
the noise-only distribution’s contour lines gives an indication 
of its significance. For instance, an actual observation with a 
percentile contour value > 99.7% would be outside 3cr of the 
noise-only distribution (under the assumption of a bivariate 
Gaussian distribution). Of our results the actual value ly¬ 
ing on the highest percentile is for J1024—0719 at 98.05% 
(or equivalently at 2.34cr from the mean of the scrambled 
data distribution). This is shown in Fig. 


16 


U nsurpr isingly, 
given the use of 100 scrambled data realisations there are 
three points further out in the distribution, and we can con¬ 
clude that this is not a significant event (although we note 
that in this case the distribution of scrambled data points 
is not quite a bivariate Gaussian as the three outlier points 
are beyond the 3(7 contour). 

The use of the SNR in these plots provides some further 
level of discrimination from interference compared to real 
signals in that data could return a high SNR signal (e.g. from 
a spectral line) that has a low odds ratio due to not matching 
the signal model or being incoherent between detectors. Such 
a situation would give an obvious outlier on plots such as 
Fig. |16| Real signals would be expected to have large values 
of both odds ratio and SNR and thus lie roughly along the 
diagonal of such plots. 


7 CONCLUSIONS 


We have investigated detection and parameter estimation 
issues for the model of gravitational wave emission from ro¬ 


tating neutron stars proposed in Jones (2010). The model is 
based on the star having a triaxial crust, coupled to an inte¬ 
rior superfluid. In the generic case of a triaxial non-aligned 
star, there is emission at both the spin frequency / and at 
2/. We have also considered two special cases, of a biaxial 
star (also with emission at / and 2/), and a triaxial aligned 
one (emission only at 2/), the last of these being the case 
conventionally assumed in gravitational wave searches. 

We have found that in the generic case of emission 
from a triaxial non-aligned star, the set of physical param¬ 


eters originally used in Jones (2010), the ‘source parame¬ 


ters’, are correlated in a highly complex way. However, a re- 
parameterisation in terms of complex waveform amplitudes 



signal-to-noise ratio 


Figure 16. The distribution of log odds ratios (for the triaxial 
non-aligned model versus noise) against SNR for 100 ‘scrambled’ 
noise-only realisations of the dataset for pulsar J1024—0719. The 
actual (i.e. not randomised in time) value is plotted as the red 
cross. Assuming that the noise-only distribution is a bivariate 
Gaussian, the 1, 2 and 3(7 (or 68.3%, 95.4% and 99.7%) probabil¬ 
ity contours are also plotted, along with the probability contour 
on which the actual value lies. 


using the ‘waveform parameters’ described in Jones (2015) 
breaks this degeneracy. When using a stochastic sampling 
method (such as nested sampling) to estimate parameter 
probability distributions from data containing such a signal, 
the complexity of the source parameter space makes a search 
there roughly half as computationally efficient as one in the 
waveform parameter space. 

For a signal described by the triaxial non-aligned model, 
we showed that estimates of many of the true individual 
source parameters, including the important parameters giv¬ 
ing the asymmetry of the moment of inertia tensor, will al¬ 
ways be poorly constrained due to the degeneracies in the 
model. This may limit the astrophysical information that 
can be extracted on such a source. We also note the (often 
overlooked) fact when discussing parameter estimation for 
these sources that for any signal there is a degeneracy in 
the full physically allowed parameter space that means the 
signal can only ever be constrained to a number of equally 
likely modes. 

Working in the waveform parameterisation, and assum¬ 
ing stationary Gaussian noise in the data, we have used sim¬ 
ulations to calculate the odds ratio for three different signal 
models compared to noise alone. We find that for a 1% false 
alarm rate, calculated from an odds ratio threshold value, 
all three models have efficiencies of close to 100% for purely 
triaxial aligned signals with SNR > 6. The simplicity of the 
triaxial aligned model compared to the others does make it 
slightly more efficient in this case, but only marginally so. 
For signals from the triaxial non-aligned and biaxial mod¬ 
els, with significant power at /, when assuming either of 
those models, efficiencies are close to 100% for SNRs > 6, 
while assuming the triaxial aligned model leads to some loss 
in efficiency, due to some proportion of signals having very 
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Table 3. Upper limits on gravitational wave amplitudes, and log odds ratios comparing the triaxial aligned model (1), biaxial model 
(2) and triaxial non-aligned model (3) signal models, for 43 isolated pulsars using LIGO S5 data. Also given is the percentile probability 
contour in the odds ratio (of triaxial non-aligned signal verses noise)—SNR plane of the background distribution on which the actual 
data point sits. 


PSR 

/ (Hz) 

2/ (Hz) 

^95% 

^21 

^95% 

/l95% 

In G>i 2 

India 

In d23 

percentile 

J0024-7204C 

173.7 

347.4 

6.5x10-26 

3.4x10-26 

6.5x10-25 

12.8 

12.8 

-0.8 

94.5 

J0024-7204D 

186.7 

373.3 

5.8x10-26 

2.0x10-26 

4.7x10-26 

13.4 

13.4 

0.7 

81.3 

J0024-7204F 

381.2 

762.3 

1.1x10-26 

4.7x10-26 

8.2x10-26 

11.1 

11.1 

0.7 

46.0 

J0024-7204G 

247.5 

495.0 

1.5x10-25 

4.2x10-26 

8.5x10-26 

12.1 

12.1 

-1.4 

59.8 

J0024-7204L 

230.1 

460.2 

6.4x10-26 

2.8x10-26 

5.9x10-26 

12.5 

12.5 

0.3 

44.7 

J0024-7204M 

272.0 

544.0 

7.2x10-26 

3.6x10-26 

7.4x10-26 

12.6 

12.6 

-0.5 

43.2 

J0024-7204N 

327.4 

654.9 

1.5x10-26 

3.2x10-26 

7.3x10-26 

12.3 

12.3 

-0.1 

58.1 

J0711-6830 

182.1 

364.2 

9.5x10-26 

2.0x10-26 

4.2x10-26 

12.1 

12.1 

-0.2 

23.2 

J1024-0719 

193.7 

387.4 

1.2x10-24 

3.3x10-26 

6.3x10-26 

11.0 

11.0 

-0.3 

98.0 

J1730-2304 

123.1 

246.2 

7.8x10-26 

3.2x10-26 

6.6x10-26 

13.1 

13.1 

1.3 

95.1 

J1744-1134 

245.4 

490.9 

6.7x10-26 

4.3x10-26 

8.3x10-26 

13.2 

13.2 

0.8 

63.8 

J1748-2446G 

118.5 

237.1 

6.2x10-26 

2.1x10-26 

4.4x10-26 

12.4 

12.4 

0.7 

43.8 

J1748-2446D 

212.1 

424.3 

3.8x10-26 

3.4x10-26 

7.5x10-26 

12.7 

12.7 

-0.3 

36.7 

J1748-2446F 

180.5 

361.0 

7.7x10-26 

3.1x10-26 

6.4x10-26 

12.2 

12.2 

0.0 

60.0 

J1748-2446H 

203.0 

406.0 

7.0x10-26 

3.5x10-26 

6.9x10-26 

12.9 

12.9 

-0.2 

61.1 

J1748-2446K 

336.7 

673.5 

2.1x10-25 

3.2x10-26 

7.0x10-26 

11.3 

11.3 

0.1 

61.6 

J1748-2446L 

445.5 

891.0 

6.3x10-26 

6.9x10-26 

1.7x10-25 

11.4 

11.4 

-0.3 

81.0 

J1748-2446R 

198.9 

397.7 

1.1x10-25 

3.8x10-26 

8.7x10-26 

12.7 

12.7 

2.2 

60.0 

J1748-2446S 

163.5 

327.0 

5.5x10-26 

2.6x10-26 

5.2x10-26 

12.6 

12.6 

-0.7 

31.5 

J1748-2446T 

141.1 

282.3 

8.3x10-26 

2.7x10-26 

5.4x10-26 

13.7 

13.7 

0.2 

54.9 

J1748-2446aa 

172.8 

345.5 

1.1x10-25 

9.9x10-26 

2.1x10-25 

12.7 

12.7 

0.1 

63.0 

J1748-2446ab 

195.3 

390.6 

6.0x10-26 

2.6x10-26 

5.1x10-26 

13.1 

13.1 

0.8 

63.1 

J1748-2446ac 

196.6 

393.2 

7.9x10-26 

2.3x10-26 

4.2x10-26 

12.8 

12.8 

-0.1 

52.9 

J1748-2446af 

302.6 

605.3 

5.5x10-25 

5.4x10-26 

1.1x10-25 

12.0 

12.0 

-3.1 

95.1 

J1748-2446ag 

224.8 

449.6 

7.5x10-26 

5.3x10-26 

1.1x10-25 

12.3 

12.3 

0.3 

90.9 

J1748-2446ah 

201.4 

402.8 

7.1x10-26 

3.1x10-26 

6.2x10-26 

12.8 

12.8 

-1.0 

74.2 

J1801-1417 

275.9 

551.7 

1.2x10-25 

2.9x10-26 

5.8x10-26 

11.7 

11.7 

-0.0 

29.4 

J1803-30 

140.8 

281.6 

6.9x10-26 

3.3x10-26 

6.5x10-26 

11.7 

11.7 

0.9 

84.5 

J1823-3021A 

183.8 

367.6 

8.7x10-26 

2.2x10-26 

4.1x10-26 

14.1 

14.1 

0.6 

82.0 

J1824-2452A 

327.4 

654.8 

8.4x10-26 

3.9x10-26 

8.4x10-26 

13.2 

13.2 

0.8 

91.2 

J1824-2452B 

152.7 

305.5 

7.7x10-26 

2.5x10-26 

5.2x10-26 

13.6 

13.6 

-0.1 

94.3 

J1824-2452E 

184.5 

369.1 

5.9x10-26 

3.5x10-26 

6.7x10-26 

12.2 

12.2 

-0.3 

47.4 

J1824-2452F 

408.0 

815.9 

1.4x10-24 

4.6x10-26 

9.1x10-26 

11.9 

11.9 

0.6 

34.5 

J1843-1113 

541.8 

1083.6 

1.8x10-25 

7.6x10-26 

1.5x10-25 

12.5 

12.5 

-0.2 

97.0 

J1905+0400 

264.2 

528.5 

6.2x10-26 

3.8x10-26 

7.8x10-26 

12.4 

12.4 

0.6 

47.8 

J1910-5959B 

119.6 

239.3 

5.9x10-26 

1.9x10-26 

4.2x10-26 

12.8 

12.8 

-0.8 

49.3 

J1910-5959G 

189.5 

379.0 

6.0x10-26 

2.2x10-26 

4.6x10-26 

12.1 

12.1 

-1.9 

70.9 

J1910-5959D 

110.7 

221.4 

7.7x10-26 

1.4x10-26 

2.9x10-26 

12.6 

12.6 

1.3 

50.6 

J1910-5959E 

218.7 

437.5 

6.3x10-26 

2.1x10-26 

4.6x10-26 

12.4 

12.4 

0.4 

32.7 

J1911+1347 

216.2 

432.3 

2.3x10-24 

2.5x10-26 

6.0x10-26 

12.2 

12.2 

-0.1 

51.7 

J1939+2134 

641.9 

1283.9 

2.4x10-25 

7.5x10-26 

1.6x10-25 

11.4 

11.4 

-0.1 

37.2 

J2124-3358 

202.8 

405.6 

8.3x10-26 

2.1x10-26 

4.7x10-26 

12.5 

12.5 

-0.4 

57.0 

J2322+2057 

208.0 

415.9 

4.3x10-26 

4.4x10-26 

8.8x10-26 

13.9 

13.9 

-0.2 

52.5 


little power at 2/. Without some idea of the true distribu¬ 
tion of signals strengths between / and 2/ we cannot say 
whether only searching a 2/ would cause any real signals to 
be missed. But, the ease of searching at both frequencies, 
and the only very minor efficiency loss, makes performing 
such a search seem sensible in the future. 

When comparing model evidences we find that for sim¬ 
ulations containing any of the three models, at very low SNR 
(< 2.5) the triaxial aligned model is always favoured. If the 
simulated signal is from the triaxial aligned model then this 
is always favoured over the other two models, whilst there is 
a strong preference for the biaxial as compared to the triaxial 
non-aligned one due to the biaxial model being the simpler 
one. For simulations containing triaxial non-aligned signals. 


the triaxial non-aligned and biaxial models are favoured over 
the triaxial aligned model for coherent SNRs > 10. We see 
a similar situation in simulations containing biaxial injec¬ 
tions: at low SNR the triaxial aligned model is favoured, but 
at higher SNR the biaxial and triaxial non-aligned models 
become favoured over the triaxial aligned one, while the bi¬ 
axial model becomes favoured over the triaxial non-aligned 
one. 

Our results show that, even though to detect all triaxial 
non-aligned (or biaxial) signals at SNR 20 one should use 
the triaxial non-aligned (or biaxial) model when computing 
the evidence, the Occam factor still significantly penalises a 
reasonable percentage of those models when deciding which 
best fits the data. As such it is worth noting that even at 
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high SNR it is often not possible to distinguish a triaxial 
non-aligned signal from a biaxial one. However, the cost of 
searching for a triaxial non-aligned signal compared to a 
biaxial signal is relatively minor, so there is no reason to 
not include such a search in the future. 

Having developed the machinery needed to search for 
such signals, we then applied our methods to real gravi¬ 
tational wave detector data. Specifically, data from the S5 
LIGO science run was used to search for two-harmonic sig¬ 
nals from 43 known pulsars with accurately known timing 
solutions, whose spin frequencies lie within the LIGO band. 
We found no gravitational wave signals, and so upper lim¬ 
its were given on the amplitude-like parameters C 21 and 
C 22 on the triaxial non-aligned model. This is the first time 
gravitational wave detector data has been searched for such 
signals. The techniques developed here will be applied to the 
more sensitive data soon to come from the new generation 
of advanced gravitational wave detectors ( Aasi et ^|2015 
Acernese et al.||2015 ) 


There is further work to do on the choice of prior proba¬ 
bility distributions that one assumes for the parameters, par¬ 
ticularly for the amplitude-like parameters. A simple choice, 
uniform up to some fixed maximum amplitude, was used 
here, but other choices are possible and will affect the results 
obtained. Glosely related to this is the issue of the physical 
interpretation of the parameters C 21 and C 22 ] it is useful to 
consider whether they have a direct physical interpretation. 
Their maximum values are presumably related to the shear 
modulus and breaking strain of the crust, and also to the 
strength of the interaction between the superfluid and non¬ 
superfluid parts of the star, but the relation is not obvious, 
and is clearly worthy of further study. 
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APPENDIX A: RELATING THE SOURCE AND 
WAVEFORM PARAMETERS 


In this appendix we briefly describe the relation between 
the source parameters and the waveform parameters, for all 
three of our chosen physical models. Full details can be found 
in Jones (2015). Note that there a few small differences in no¬ 


tation between the equations given here and those of | Jones] 


(2015). The angle -0 of 


Jones 


(2015) is denoted by A in this 


paper, as we reserve the symbol '0 for the polarisation an¬ 
gle, describing the projection of the star’s spin axis on the 
plane of the sky. Also, the quantities denoted by {€ 22 ^ 021 ) 
of Jones (2015) are denoted by (C 22 ,C 2 i) in this paper. 
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A1 Triaxial non-aligned case 

For the triaxial non-aligned case, the waveform is given in 
terms of source parameters by Eqns 0 and 0- and by 
Eqns 0 and 0 in terms of waveform parameters. It can 
be shown that the relationship between the five source pa¬ 
rameters (/ 21 ,/ 31 , 00 , A) and the four source parameters 

(^ 21 , ^ 21 , ^^ 22 , ^ 22 ) is as follows: 

C 22 = 2 {[/ 2 i(sin^ A — cos^ Acos^ 0 ) — I 31 sin^ 0 ]“^-\- 

(/21 sin 2A cos (Al) 

^^21 = 2 {(/ 2 i sin 2A sin 0)^ + (/21 cos^ A — /ai)^ sin^ 20}^^"^, 

(A 2 ) 


tan[ 20 o - ^ 22 ] 


tan [00 - ^ 21 ] 


I 21 sin 2 Acos^ 
I 21 (sin^ A — cos^ A cos^ 0) 

(/21 cos^ A — / 31 ) sin 20 
/2isin2Asin^ ’ 


— /31 sin^ 6 ’ 

(A 3 ) 

(A 4 ) 


see equations (62)-(65) of 


Jones 


(2015). Given the set of 


five source parameters one can calculate unique values for 
the four waveform parameters. However, for a set of four 
waveform parameters, the corresponding solution for the five 
source parameters will have one degree of freedom, generat¬ 
ing the sorts of complex structure in source parameter space 
seen in Eigs[^and[^ 


A2 Biaxial case 

If we set I 21 = 0 in Eqns 0 and 0 we obtain the biaxial 


signal in terms of source parameters: 

— — 2 / 31(1 + cos^ i) sin^ 0 cos 2(0.1 + 0 o), (A5) 

= — 4/31 cos ^sin^ ^ sin 2(Qt + 0o), (A6) 

— — /31 sin L cos l sin 20 sm(Ot + 0 o), 

h-x =^31 sin^sin2^cos(Qt + 0o), (A 8 ) 


where we have separated out the ‘+’ and ‘x’ polarisation 
components. 

It can be shown that the corresponding waveform pa- 
rameterisation can be written as Eqn. 0 and ([^, with the 
extra condition 


^22 — 24 >? i : 


(A 9 ) 


see equation (90) of Jones (2015). The relation between the 
source and waveform parameters can be shown to be 


iW) 


hi = — 2^22 
00 


1 + 


G 21 V 


see equations (82), ( 88 ) and (89) of 


Jones 


(2015). 


(AlO) 

(All) 

(A 12 ) 


A3 Triaxial aligned case 

If we set ^ = 0 in Eqns 0 and 0 we obtain the triaxial 
aligned signal in terms of source parameters. There in no 


signal at frequency /, leaving only: 

= — 2 / 21(1 + cos^ ^) cos2[Qt + (00 + A)], (A13) 

= — 4/21 cos ^ sin 2 [fit + (00 + A)]. (^14) 


The angles 0o and A are now degenerate, with only their 
sum appearing in the waveform. This sum can now be taken 
as replacing the separate parameters, or else one can simply 
set one of them to a constant value and search over the other 
(e.g. set A = 0 and search over 0o). 

The corresponding waveform parameterisation has 
C 21 = 0 so that we only have the 2/ signal given by Eqn. 0 . 
with the source and waveform parameters being related by 

hi =^022, (A15) 

2(00 + A) (A16) 


see equations (103) and (104) of |Jones| ( [m^ . 


APPENDIX B: ALGORITHM TO FILL IN THE 
FULL SOURCE PARAMETER SPACE 

In this appendix we present (as pseudo-code) the algorithm 
used to map the minimal source parameters space, whose 
ranges were given in Table to the full parameter space. 


Algorithm 1 Transformations required to map from the 
minimal range in the source model parameter space to the 
full space. The values with the superscript ‘min’ represent 
the original samples from the minimal parameter range. 


1 : 

2 : 


3 

4 

5 

6 

7 

8 
9 

10 


for i ^ 1 to 4 do 

draw a random unique set consisting of one quarter 
of the samples, such that j- is a vector of those samples’ 
indices for loop i 

SUi) = + {i- l)7r/2 

0Ui) = 

MJi) = 

HJi) = 
hlij,) = 
hiUi) = lEiGJi) 

for k ^ 1 to i do 


Onew = arctan 


sine(ji) sin A(jf^) 


11 : MJi) = 

f cos(/)o(jf^)cosA(jf^)-sin(^o(jf^)cos0(jf^)sinA(jf^) \ 

- cos cfo (j cos e (jf i ) sin A (jf i ) - sin cfo (j ^ ) cos A (j ^ ) y 

12: \{ji) = arctan { ,infu°)!osX(jo ) 

13: hi{ji) = hl{j^] - hi(ji) 

14: hi(ji) = hi{ji) 

15: Oi — ^new 

16: end for 

17: end for 

18: randomly select half the samples, such that fci is a vector 
of those sample indices 
19: apply X(ki) = X(ki) + tt 

20: randomly select half the samples, such that k 2 is a vector 
of those sample indices 

21: apply 0 (k 2 ) = tt — ^(^ 2 ), 00 (^ 2 ) = 00 (^ 2 ) + tt and 

\(k2) = -X(k2) 

22 : 00 = 00 (mod 2n) 

23: A = A (mod 2ti) 

































